/*****************************************************************************/
/*                                                                           */
/*    Domino Operation System Mathematics Module                             */
/*                                                                           */
/*    Copyright (C) 2005 Laszlo Arvai                                        */
/*                                                                           */
/*    ------------------------------------------------------------------     */
/*    mathFFT.h - Fast fourier transformation routines                       */
/*****************************************************************************/

#include "mathComplex.h"
#include "mathFFT.h"

///////////////////////////////////////////////////////////////////////////////
// Module global variables
static dosComplex l_fft_buffer[FFT_NUMBER_OF_POINTS/2];

///////////////////////////////////////////////////////////////////////////////
// FFT routine
void mathFFT( dosInt32* in_data )
{
	dosInt16 p,p2,q,r,i,nper2p,nper4p;
  dosInt32 s,*poi1,*poi2;
  dosComplex *p_exp;
  dosComplex *p_from,*p_to,*p1_from,*p1_to,*p2_from,*p2_to;
  dosComplex cx,cxexp;

  /***********************   FIRST CYCLE   ********************/

  p = FFT_NUMBER_OF_POINTS/2;
  poi1 = in_data;
  poi2 = in_data+p;
  p1_to = l_fft_buffer;
  for(q = 0; q < p; q++)
  {
		p1_to->re = *poi1+(*poi2);
    p1_to->im = *poi1-(*poi2);
    poi1++;
    poi2++;
    p1_to++;
  }

  /******************  MAIN LOOP BEGINNING   ******************/

  p_from = l_fft_buffer;
  p_to = (dosComplex *)in_data;
  for(i = 1; i < FFT_LOG_NUMBER_OF_POINTS; i++)
  { 
		p2 = p;
    p >>= 1;
    nper2p = FFT_NUMBER_OF_POINTS / p2;

    /****************   CASE r = 0 & r = FFT_NUMBER_OF_POINTS/2p   *******************/

    p2_from = p+(p1_from = p_from);
    p1_to = p_to;
    for(q = 0;q<p;q++) { p1_to->re = p1_from->re+p2_from->re;
		       p1_to->im = p1_from->re-p2_from->re;
		       p1_from++;
		       p2_from++;
		       p1_to++;
		     }
    p1_from += p;
    p2_from += p;
    p2_to = p_to+FFT_NUMBER_OF_POINTS/2-p;

    /*****************   GENERAL  CASE   **********************/

    if( FFT_NUMBER_OF_POINTS >= 4*p )
		{
			p_exp = g_fft_exp_table + p - 1;
			nper4p = nper2p/2;
			for(r = 1; r < nper4p; r++)
	    { 
				cxexp = *p_exp;
				for(q = 0;q<p;q++)
				{
					mathCplxMultFFT( p2_from, &cxexp, &cx);
					mathCplxAdd(p1_from,&cx,p1_to);
					mathCplxSub(p1_from,&cx,p2_to);
					p2_to->im = -p2_to->im;
					p1_from++;
					p2_from++;
					p1_to++;
					p2_to++;
				}
	      p_exp += p;
	      p1_from += p;
	      p2_from += p;
	      p2_to -= p2;
	    }
     }

    /********************   r = FFT_NUMBER_OF_POINTS/4p case   *********************/

    p2_from = p+(p1_from = p_from);
    for(q = 0; q < p; q++)
    {
			p1_to->re = p1_from->im;
			p1_to->im = -p2_from->im;
			p1_from++;
			p2_from++;
			p1_to++;
    }

    /***************   END OF THE MAIN LOOP   *****************/
    p1_from = p_from;
    p_from = p_to;
    p_to = p1_from;
  }

  if(FFT_LOG_NUMBER_OF_POINTS%2)
  {
		poi1 = (dosInt32*)l_fft_buffer;
    poi2 = in_data;
    for(i = 0;i<FFT_NUMBER_OF_POINTS;i++)
		  *(poi2++) = *(poi1++);
  }

  s = (dosInt32)FFT_NUMBER_OF_POINTS;
	for(i = 0;i<FFT_NUMBER_OF_POINTS;i++) 
		*(in_data++) /= s;
}

